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1  SUMMARY 

Ground-based  optical  and  radar  sites  routinely  acquire  resolved  images  of  satellites.  These  images  provide  the 
means  to  construct  accurate  wire-frame  models  of  the  observed  body,  as  well  as  an  understanding  of  its  orientation 
as  a  function  of  time.  Unfortunately,  because  such  images  are  typically  acquired  at  a  single  wavelength,  this  kind  of 
analysis  provides  little  or  no  information  on  the  types  of  materials  covering  the  satellite’s  various  surfaces.  Detailed 
surface  material  characterization  generally  requires  multi-band  radiometric  and/or  spectrometric  measurements. 
Many  widely-available  instruments  provide  such  multi -band  information  (e.g.,  spectrographs  and  multi-channel 
photometers).  However,  these  sensors  typically  measure  the  brightness  of  sunlight  reflected  from  the  entire  satellite, 
with  no  spatial  resolution  at  all.  Because  such  whole -body  measurements  represent  a  summation  of  contributions 
from  many  reflecting  surfaces,  an  “un-mixing”  analysis  must  be  employed  to  characterize  the  reflectance  of  the 
satellite’s  individual  sub -components.  The  objective  of  this  paper  is  to  outline  the  theory  required  for  such  an  un¬ 
mixing  process,  focusing  on  two  newly-developed  analysis  methods.  Both  methods  retrieve  satellite  surface 
properties  from  temporal  sequences  of  whole -body  brightness  measurements.  Both  require  the  following  as  input: 
1)  a  set  of  multi -band  measurements  of  a  satellite’s  brightness  in  reflected  sunlight,  2)  the  satellite’s  wire-frame 
shape  model,  including  each  major  exterior  sub-component  capable  of  reflecting  sunlight,  3)  the  satellite’s  attitude, 
specifying  the  orientation  of  all  of  the  body’s  components  at  the  times  of  each  measurement.  The  first  method  also 
requires  a  library  of  bi-directional  reflection  distribution  functions  (BRDFs)  for  a  set  of  candidate  materials  covering 
the  satellite’s  surfaces,  and  yields  estimates  of  the  fraction  of  each  satellite  sub-component  covered  by  each 
candidate  material.  The  second  method  does  not  require  pre-tabulated  BRDFs,  but  instead  attempts  to  retrieve 
BRDFs  for  each  satellite  sub-component  from  the  non-resolved  data  using  a  series  expansion  approach. 

This  research  funded  by  the  Air  Force  Office  of  Scientific  Research. 

2  INTRODUCTION 

Ground-based  optical  and  radar  sites  routinely  acquire  resolved  images  of  satellites,  yielding  a  great  deal  of 
knowledge  about  orbiting  spacecraft.  In  particular,  the  Air  Force  Maui  Optical  and  Supercomputing  (AMOS) 
Detachment  on  Maui  has  been  acquiring  optical  imagery  using  two  work-horse  imagery  systems.  The  AMOS 
Advanced  Electro-optical  System  (AEOS)  3.6m  telescope  provides  visible-band  and  long- wavelength  thermal 
infrared  images  with  adaptive  optics  compensation  to  remove  atmospheric  blurring.  In  addition,  the  AMOS  Gemini 
1.6m  telescope  system  provides  daytime  visible -band  and  near-infrared  speckle  images.  These  systems  reveal  a 
great  deal  of  resolvable  detail  for  satellites  in  low-Earth  orbit  (LEO),  especially  after  the  data  undergo  post¬ 
processing  enhancement  at  the  AFRL  Maui  High  Performance  Computing  Center.  Other  observatories,  such  as 
Starfire  Optical  Range  (SOR),  as  well  as  some  ground-based  radar  sites,  also  acquire  images  of  comparable  quality. 
From  such  images,  detailed  wire-frame  models  of  the  observed  satellites  can  be  assembled  and  aligned  to  the 
images.  This  process  essentially  translates  the  two-dimensional  imagery  into  detailed  three-dimensional  information 
about  the  sizes,  shapes,  and  relative  orientations  of  various  spacecraft  components.  Unfortunately,  such  image 
analysis  procedures  provide  little  or  no  information  on  the  material  and/or  optical  properties  of  the  satellite  surfaces, 
because  the  wire-frame  models  are  typically  based  on  single  spectral  band  imaging  data.  Detailed  surface  material 
and  property  characterization  generally  requires  multi-band  photometric  and/or  spectrometric  measurements. 


3  ANALYSIS  FORMULATION 

The  objective  of  this  paper  is  to  outline  the  theory  required  to  retrieve  satellite  surface  properties  from  temporal 
sequences  of  whole-body,  multi-band  brightness  measurements,  focusing  on  two  newly-developed  analysis 
methods.  The  first  method,  referred  to  here  as  “material  abundance  estimation”  (MAE),  seeks  to  determine  the 
fractional  abundances  of  materials  covering  the  major  exterior  components  of  the  observed  satellite.  This  method 
has  previously  been  described  in  detail  and  applied  to  simulated  data  of  a  convex  object  [1].  The  second  method, 
referred  to  here  as  “reflectance  distribution  estimation”  (RDE),  seeks  to  determine  the  bi-directional  reflectance 
distribution  functions  (BRDFs)  of  the  satellite’s  components. 
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Both  analysis  methods  require  three  main  types  of  input  information:  1)  a  set  of  multi -band  measurements  of  a 
satellite’s  reflected- sunlight  brightness,  2)  the  satellite’s  wire-frame  model,  including  each  major  exterior  sub¬ 
component  capable  of  reflecting  sunlight,  3)  the  satellite’s  attitude,  specifying  the  orientation  of  all  of  the  body’s 
components  at  the  time  of  each  multi-band  measurement.  The  second  two  of  these  required  inputs  can  be  derived 
through  a  variety  of  means,  including  open-source  publications  or  analysis  of  resolved  imagery  of  the  satellite 
(either  pre-launch  or  on-orbit),  as  described  previously  [1].  The  first  input  constitutes  the  data  to  be  analyzed,  and 
comprise  observations  from  one  (or  more)  sensors  that  provide  non-resolved,  multi-band  brightnesses  of  the 
satellite.  Such  measurements  naturally  contain  an  abundance  of  information  on  the  material -specific  reflective 
properties  of  the  satellite’s  surfaces.  However,  such  spectrometric  instruments  often  provide  brightness 
measurements  of  the  entire  object ,  rather  than  from  individual  satellite  surfaces  or  component.  In  other  words,  each 
measurement  represents  a  summation  of  light  reflected  from  many  satellite  components,  thereby  requiring  an  “un¬ 
mixing”  analysis  process  to  separate  and  retrieve  individual  component  reflectance  properties. 

Because  the  first,  MAE  method  seeks  to  estimate  material  abundances,  it  requires  one  additional  input:  pre-tabulated 
BRDFs  for  a  set  of  candidate  materials  covering  the  satellite’s  surfaces.  It  then  yields  as  output  estimates  of  the 
fraction  of  each  major  satellite  sub-component  covered  by  each  material.  The  MAE  method  can  suffer  significantly, 
however,  when  provided  with  a  database  not  containing  materials  actually  covering  one  or  more  components  of  the 
observed  satellite  [1],  or  not  containing  information  on  how  the  BRDFs  change  as  materials  experience  space 
weathering  effects.  These  two  BRDF  database  limitations  can  potentially  significantly  limit  the  effectiveness  of  the 
MAE  method  when  applied  to  unknown  or  aging  satellites.  The  second,  RDE  method  was  developed  specifically  to 
address  these  limitations.  It  does  not  require  any  material  BRDF  information  as  input.  Instead,  the  RDE  method 
attempts  to  retrieve  BRDF  parameters  for  each  satellite  sub-component  directly  from  the  non-resolved  data.  For  this 
reason,  the  RDE  method  promises  to  be  more  appropriate  for  analyses  of  objects  covered  with  unknown  or  exotic 
materials  or  whose  reflectance  has  changed  significantly  under  the  effects  of  space  weathering. 

3.1  Satellite  Wireframe  Shape  Models 

The  required  satellite  wireframe  models  can  be  assembled  by  combining  one  or  more  data  sources:  design 
information  provided  by  the  manufacturer,  pre-launch  photographs  and/or  artistic  renderings,  and/or  from  resolved 
imagery  of  already-orbiting  objects  [1].  Wireframe  models  are  typically  assembled  using  primitive  components 
(such  as  flat  panels,  spheres,  cylinders,  cones,  parabolic  dishes,  etc.)  combined  to  form  accurate  three-dimensional 
representations  of  all  major  exterior  components  of  the  satellite.  After  such  model  development  is  complete,  each 
component  of  the  wire-frame  model  can  then  be  decomposed  into  a  series  of  perfectly  flat  facets.  Even  round 
components  may  be  approximated  in  this  manner  using  many  small  facets.  Each  facet  comprises  a  planar  polygon 
with  shape,  area,  and  orientation  derived  from  the  original  wire -frame  model. 

3.2  Satellite  Attitude  Models 

Satellite  attitude  model  can  also  be  assembled  in  several  ways.  For  instance,  cooperative  satellite  owner/operators 
may  willingly  provide  “quaternion  data”  encoding  the  detailed  body  attitude  as  a  function  of  time,  or  at  least  provide 
the  attitude  control  system  (ACS)  stabilization  parameters  employed  on-board  the  satellite  during  the  period  of 
observations.  Alternatively,  analysis  of  high-quality  images  can  provide  the  orientation  of  on-orbit  satellites  as  a 
function  of  time,  derived  from  the  sequential  frame-to-frame  adjustments  required  to  align  the  wire-frame  model  to  a 
series  of  observed  images  [1].  Detailed  image  analysis  can  reveal  ACS  operational  modes,  or,  alternatively,  indicate 
rotational  parameters  for  spin- stabilized  satellites  or  even  tumbling  objects.  Ultimately,  this  information  can 
conceivably  be  used  to  assemble  a  predictive  attitude  model,  and  used  as  input  by  both  of  the  analysis  methods 
discussed  here  to  specify  the  body’s  orientation  at  the  time  of  each  multi-band  measurement. 

In  mathematical  terms,  the  attitude  model  comprises  the  set  of  parameters  required  to  specify  the  body’s  “attitude 
matrix,”  R,  as  a  function  of  time.  This  3x3  rotation  matrix  converts  vectors  from  an  inertial  reference  frame  (taken 
here  to  be  the  Earth-centered  J2000  frame)  into  the  body-fixed  reference  frame,  written  symbolically  as  follows: 

R  =  R(0  (i) 

where  the  time  dependence  indicates  the  changing  orientation  of  the  satellite.  As  an  example,  the  attitude  matrix 
provides  the  means  to  convert  the  inertial-frame  satellite-to-observer  and  unit  direction  vector,  o  (t),  and  the 
satellite-to-Sun  unit  vector,  s  (t),  into  the  body-frame  through  matrix  multiplication: 

o(t)  =  R(t)  o \t)  s(/)  =  R(/)  s *(t)  (2) 
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In  this  discussion,  the  superscript  denotes  inertial-frame  vectors,  whereas  body-frame  vectors  are  denoted 
without  superscripts.  The  two  vectors  given  in  Eq.  (2)  play  an  important  role  in  determining  the  observed  brightness 
of  sunlight  reflected  from  the  body,  as  described  below.  While  both  of  these  vectors  depend  on  time,  in  this 
discussion  their  explicit  time  dependence  will  often  be  suppressed  for  brevity. 

3.3  Theory  for  Multi-band  Un-mixing  Analysis 

The  satellite  wire-frame  and  attitude  models  described  above  provide  a  means  of  calculating  the  visibility  and 
observation  geometry  of  each  resolved  spacecraft  surface  during  observations  from  any  sensor  with  a  known 
location.  Notably,  if  one  also  had  detailed  a  priori  knowledge  of  the  optical  properties  of  each  surface  —  i.e.,  the 
material  composition,  and  associated  material-specific  bi-directional  reflectance  distribution  functions  (BRDFs)  — 
then  a  forward  model  of  the  spectral  signature  of  the  entire  body  could  be  constructed  by  summing  contributions 
from  all  reflecting  surfaces  visible  to  the  sensor.  However,  in  the  absence  of  such  a  priori  knowledge, 
characterizing  the  surface  properties  requires  an  un-mixing  analysis  of  the  whole-body  multi-band  measurements. 
Ideally,  these  measurements  would  be  obtained  by  spectrometric  instrumentation  with  good  spectral  resolution  and  a 
quick  cadence.  This  would  provide  the  spectral  data  required  to  discriminate  different  surface  materials,  and  the 
temporal  data  to  distinguish  shiny  and  dull  surfaces.  The  same  task  could  also  conceivably  be  performed  by  multi¬ 
channel  photometer(s)  providing  relatively  broad-band  spectral  coverage,  albeit  less  efficiently. 
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Figure  1.  BRDF  plots  for  three  common  spacecraft  materials  tabulated  in  the  TASAT  software  package  [2].  The  top  panels  plot  BRDF  surface 
representations  at  a  wavelength  of  X  =  0.6  pm,  illustrating  relative  contributions  of  specular  and  diffuse  reflection  components.  The  yellow  line 
indicates  an  example  illumination  direction,  incident  onto  a  reflecting  sample  of  material  lying  in  the  bottom,  horizontal  plane.  The  green  line 
shows  the  corresponding  direction  of  specular  reflection.  Specular  (or  shiny)  reflections  appear  as  long  and  sharp  spikes  along  this  green  line  in 
these  surface  plots.  The  bottom  panels  plot  the  total  hemispherical  reflectance  (i.e.,  albedo)  spectra,  spanning  wavelengths  0.4  <  X  <  1.2  pm. 


3.3.1  The  Spectral  Intensity  of  a  Satellite  Reflecting  Sunlight 

As  discussed  earlier,  the  wire-frame  models  provide  the  orientation  and  area  of  each  facet  of  each  satellite 
component.  Specifically,  the  kth  facet  of  the  /h  component  may  be  characterized  by  its  surface  area,  Aj>h  and  normal 
unit  vector,  nj>k.  These  facets  could,  in  principle,  have  time  dependent  orientations,  such  as  those  on 

articulating  spacecraft  structures,  but  this  potential  time  dependence  will  be  suppressed  here  for  brevity.  The 
spectral  intensity  (or  spectral  radiance),  L,  of  sunlight  reflected  from  the  entire  object  has  units  of  W  ster'1  pm'1.  It  is 
a  function  of  both  time  and  wavelength,  and  may  be  expressed  as  a  summation  over  the  component  facets  [1]: 


L(t,A)  =  Fs(t,A.) 


o  n 


S )  pjO l,nt,o,s)  yM(o,  s) 


(3) 


where  FSun(t,  1)  denotes  the  illuminating  solar  irradiance  (W  m‘2  pm _1),  o  and  s  denote  the  time-dependent  satellite- 
to-observer  and  satellite-to-Sun  unit  vectors  from  Eq.  (2),  and  the  function  pj  (2,  n,  o,  s)  denotes  the  surface  BRDF 
for  satellite  component  j.  Angular  brackets  denote  the  non-negative  operator,  <v>=max(0,v),  and  the  non-negative 
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dot  products  ensure  that  contributions  arise  only  from  facets  showing  an  illuminated  side  to  the  observer.  The 
shadowing/obscuration  (SO)  function,  T^o,  s),  denotes  the  fraction  of  each  facet  that  is  not  shadowed  nor  obscured 
by  other  satellite  surfaces.  For  convex  bodies  =  1  for  all  facets.  For  non-convex  bodies  this  function  generally 
varies  with  time.  The  relative  positions  and  sizes  of  the  satellite  components  specified  in  the  wire-frame  model, 
combined  with  the  o(t)  and  s (t)  unit  vectors,  provide  a  means  to  calculate  the  SO  function  using  ray-tracing,  z- 
buffering,  or  similar  algorithms  (the  currently-implemented  MAE  and  RDE  methods  employ  z-buffering). 

3.3.2  Material  Bi-directional  Reflectance  Distribution  Functions  (BRDFs) 

Notably,  the  satellite  wire-frame  and  attitude  models  together  provide  all  of  the  quantities  required  to  calculate  the 
spectral  intensity  using  Eq.  (3),  except  for  the  BRDFs.  In  fact,  the  importance  of  the  wavelength-dependent  BRDF 
in  Eq.  (3)  cannot  be  overemphasized  in  this  regard.  It  basically  indicates  how  multi-band  measurements  can  be  used 
to  diagnose  the  material  composition,  and  is  fundamental  to  the  feasibility  of  the  un-mixing  analyses.  Many  groups 
have  measured  material  BRDFs  in  laboratory  environments  and/or  created  numerical  BRDF  models.  For  instance, 
BRDFs  for  several  spacecraft  materials  —  such  as  solar  array  panels,  milled  aluminum,  anodized  aluminum,  multi¬ 
layer  insulation,  white  paint,  etc.  —  are  available  as  part  of  the  TASAT  satellite  radiometry  simulation  software 
package  [2].  Fig.  1  illustrates  BRDFs  for  three  common  satellite  materials  tabulated  in  the  TASAT  database  [1]. 
Material  BRDFs  contain  two  general  types  of  information  that  can  help  identify  and  discriminate  satellite  materials. 
First  BRDFs  contain  spectral  “fingerprints”  (Fig.  1,  bottom  panels).  These  are  wavelength-dependent  patterns  in  the 
reflectance  (i.e.  albedo)  that  can  uniquely  identify  the  material  composition  [3,  4].  BRDFs  for  different  materials 
also  vary  significantly  in  their  relative  fractions  of  specular  vs.  diffuse  (i.e.,  shiny  vs.  dull)  reflectance.  The  two 
analysis  methods  presented  here  exploit  both  of  these  aspects  of  BRDFs.  The  MAE  method  focuses  using  a 
database  of  candidate  material  BRDFs  to  estimate  the  abundances  of  materials  on  exterior  satellite  components.  The 
RDE  method  seeks  to  estimate  the  BRDF  of  each  satellite  component  using  a  series  expansion  approach. 

3.3.3  Analysis  Method  1:  Material  Abundance  Estimation  (MAE) 

As  a  first  step,  the  MAE  formulation  assumes  that  the  surfaces  of  each  satellite  component  can  be  modeled  as  a 
mixed  set  of  distinct  materials  compiled  in  a  BRDF  database  [1].  In  this  case,  the  effective  BRDF  for  the  kth  facet  of 
the  jth  satellite  component  may  be  written  as  a  sum  over  the  individual  material  BRDFs  as  follows: 

Pj  (A,  njk ,  o,  s)  =  fUm  pm  (A,  njk ,  o,  s )  (4) 


where  fj>m  denotes  the  fractional  area  of  component  j  covered  by  material  m,  and  (1,  n,  o,  s)  denotes  the  pre¬ 
tabulated,  laboratory-measured  BRDF  for  the  pure  material  m.  Inserting  Eq.  (4)  into  Eq.  (3)  and  re-arranging  yields: 

=  (5) 

j,m 


where  the  kernel  function  is 


Kjm(t,A)  =  FSun(f„W 2 A*  (nM  -°)  (nM  -s)  nM’0’  s) s) 


(6) 


Notably,  this  kernel  can  be  readily  calculated  because  all  of  its  component  quantities  are  known  from  either 
independent  measurement  (such  as  the  solar  flux  and  the  database  of  laboratory-measured  BRDFs)  or  from  the  wire¬ 
frame  and  attitude  models.  The  only  remaining  unknown  quantities  in  Eq.  (5)  are  the  fractional  areas  (i.e., 
“abundances”)  covered  by  the  pure  materials,  fj>m.  These  quantities  are  subject  to  the  following  two  constraints: 

o<4,„<i  £/,„  =  1  (7) 

m 


The  objective  of  the  MAE  process  is  to  find  the  set  of  fractional  areas,  fj>m,  that  best  reproduce  the  observed  spectral 
radiance  data  but  that  also  satisfy  these  two  constraints. 

Typical  long- slit  spectrographs  or  broad-band  photometric  instruments  do  not  provide  continuous  measurements  of 
the  spectral  intensity,  Lit ,  X),  but  instead  provide  observations  at  a  discrete  set  of  times,  {th  i=\  ...A/}  and 
wavelengths,  {kh  1=1...  N]}.  This  formulation  idealizes  each  of  these  measurements  as  instantaneous  (i.e., 
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neglecting  the  finite  exposure  time  spanned  by  each  measurement)  with  perfectly  narrow  wavelength  sampling  (i.e., 
neglecting  the  finite  width  of  each  spectral  channel).  With  these  assumptions,  the  spectrometric  measurements  can 
be  organized  into  a  discrete  matrix  as  follows: 

L;  /  =L(/,a/)  (8) 


and  Eq.  (5)  can  be  used  to  write  the  system  of  equations  that  must  be  solved  in  the  inversion  process 

^ i,l  ^  \f  j,m  ^ i,l,j,m 

j,m 


(9) 


At  this  point,  it  is  convenient  to  combine  indices  to  streamline  the  formulation.  The  two  observation  indices  (/,  /) 
can  be  combined  into  one  master  index,  p ,  spanning  p-  1 . .  .Nj  Nj.  Similarly,  the  indices  (j,m)  can  be  combined  into 
a  master  index  v.  Using  these  master  indices,  Eq.  (9)  can  be  re-written  in  a  relatively  simple  form: 


■'Lf- 


K 


(10) 


At  first  glance,  Eq.  (10)  appears  to  be  a  linear  system,  which  could  be  solved  using  a  variety  of  methods  such  as  a 
least-squares  analysis  using  singular- value  decomposition  or  using  the  pseudo-inverse  matrix  formalism  [5]. 
However,  the  constraints  stated  in  Eq.  (7)  still  must  be  satisfied,  making  such  methods  inappropriate.  Fortunately, 
there  are  efficient  numerical  methods  [6,  7]  that  do  provide  a  means  of  imposing  such  constraints  (pre-programmed 
into  the  Matlab  software  system  as  function  Isqlin).  Solving  Eq.  (10)  using  such  methods  provides  the  set  of  fv  or, 
alternatively,  fj>m,  which  represent  the  best-fit  fractional  areas  of  each  satellite  component  covered  by  each  reflective 
material. 

3.3.4  Analysis  Method  2:  Reflectance  Distribution  Estimation  (RDE) 

The  RDE  method  does  not  assume  that  the  satellite  components  can  be  modeled  as  a  mixture  of  known  materials,  as 
does  the  MDE  method.  Instead,  the  RDE  method  attempts  to  estimate  the  BRDF  of  each  satellite  component  using  a 
series  expansion  approach.  In  this  case,  the  model  BRDF  for  the  kth  facet  of  the  /h  satellite  component  may  be 
written  as  a  sum  over  a  set  of  BRDF  expansion  functions  as  follows: 

PjU,  njk ,  o,  s)  =  Yuaj,m  (ll) 


where  aj>m  denotes  the  albedo  of  component  j  associated  with  the  BRDF  expansion  function  (2,  n,  o,  s).  These 
BRDF  expansion  functions  do  not  correspond  to  pure  materials,  but  instead  are  analytical  functions  that  span  a  range 
of  reflectance  distributions  from  very  dull  (diffuse)  to  very  shiny  (specular),  as  illustrated  in  Fig.  2.  In  this 
formulation,  these  expansion  functions  have  the  following  form: 


/U  n,o,s)  = 


J  1/7T 
\/hrx  (n,  O,  S,  ) 


m  =  0 

m  —  1  •  •  *^SPEC 


(12) 


The  m  =  0  component  represents  purely  diffuse  reflection  using  the  simple  Lambertian  BRDF.  The  remaining 
components  represent  increasingly  specular  reflection  using  analytical  functions  based  on  the  “Cook-Torrance” 
formulation  [8].  This  series  expansion  approach  assumes  that  real-world  BRDFs,  such  as  the  surfaces  shown  in  Fig. 
1,  can  be  approximated  sufficiently  using  a  weighted  sum  of  analytical  functions,  like  those  in  Fig.  2. 

One  choice  for  series  expansion  functions  is  based  on  the  Cook-Torrance  BRDF  formulation  [8],  used  here  because 
it  is  both  physics-based  and  relatively  computationally  efficient.  The  underlying  assumption  for  this  formulation  is 
that  most  surfaces  are  uneven  at  a  size  scale  comparable  to  or  larger  than  the  wavelength  of  reflected  light.  BRDFs 
for  such  “rough”  surfaces  must  account  for  the  collective  reflection,  shadowing  and  obscuration  from  the  “micro¬ 
facet  distributions”  of  the  surfaces  [9,  10],  and  have  effective  specular  components  that  form  broad  distributions 
peaked  in  the  direction  of  perfect  mirror-like  reflection.  The  exact  form  of  these  broadened  peaks  depends  on  the 
parameters  of  the  micro-facet  distributions.  However,  the  Cook-Torrance  formulation  [8]  uses  a  single,  physics- 
based  parameter  —  p,  the  RMS  slope  of  the  micro-facets  —  to  specify  the  width  of  the  specular  peak  as  follows: 
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Pc'x  (n,  o,  s,  fJ) 


1 

;r(n-o)(n-s) 


^-(tan  at  juf 

4//  cos 4  a 


f 


min  1, 


v 


2(n  b)(n  o) 

(ob) 


2(n*b)(n*s)^ 
(o-b)  ; 


(13) 


where  a  the  angle  between  the  surface  normal,  n,  and  the  vector  b,  which  bisects  the  o  and  s  vectors.  The  first  two 
factors  in  Eq.  (13)  account  for  the  collective  reflection  from  the  micro-facet  distribution  and  the  last  for  micro-facet 
shadowing/obscuration  effects  [8],  not  to  be  confused  with  the  macroscopic  facet  SO  function  introduced  in  Eq.  (3). 
The  set  {pm}  in  Eq.  (12)  should  span  the  range  of  RMS  slopes  representative  of  common  materials;  something  like 
0.04  <fim<  0.5  [see  8]. 


Figure  2.  Surface  plots  for  six  BRDF  expansion  functions  used  in  the  RDE  method,  spanning  diffuse  reflection  (upper  left)  to  strongly  specular 
reflection  (lower  right).  In  each  plot,  the  yellow  line  indicates  an  incident  illumination  direction,  and  the  green  line  the  corresponding  direction  of 
specular  (i.e.,  mirror-like)  reflection.  Strongly  specular  BRDFs  are  dominated  with  long,  sharp  spikes  along  this  green  line  (lower  right  panel). 
Diffuse  BRDFs  don’t  have  such  directionality;  in  fact,  perfect  Lambertian  reflectors  have  hemi-spherical  surfaces  (upper  left). 

Inserting  Eq.  (11)  into  Eq.  (3)  and  re-arranging  yields: 

L(t,X)  =  Krj  m(t,X)  (i4) 

j,m 


where  this  kernel  function  has  similar  form  to  that  given  in  Eq.  (6),  and  also  can  be  readily  calculated  because  all  of 
its  component  quantities  are  known  from  either  independent  measurement,  the  wire-frame  and  attitude  models,  or 
from  the  BRDF  expansion  functions.  The  only  remaining  unknown  quantities  in  Eq.  (14)  is  the  set  {aym},  which 
denote  the  albedo  of  component  j  associated  with  the  BRDF  expansion  function  m.  These  are  subject  to  the 
following  two  constraints: 


a  =/  a- 

m 


<  1 


(15) 


The  objective  of  the  RDE  analysis  is  to  find  the  set  of  albedos  {ajjm}  that  best  reproduce  the  spectral  radiance  data 
but  that  also  satisfy  these  two  constraints.  Note  that  the  second  of  these  constraints  differs  from  the  MAE  method, 
in  that  it  is  summation  equality  constraint  rather  than  a  summation  inequality  constraint.  Fortunately,  both  can  be 
handled  by  available  numerical  methods  [6,  7].  Solving  for  the  best- fit  set  {aj>m}  can  formulated  and  performed  in  a 
fashion  very  similar  to  the  MAE  method  described  above.  Eqs.  (11)-(13)  can  be  then  used  to  calculate  the  best 
estimate  of  the  BRDF  for  each  satellite  component,  pj  {X,  n,  o,  s). 

4  APPLICATION  TO  SIMULATED  OBSERVATIONS 

This  section  applies  the  two  analysis  method  to  simulated  observations  of  a  relatively  simple  non-convex  model 
satellite,  as  shown  Fig.  3  and  detailed  in  Table  1.  Simulated  spectral  irradiances  were  calculated  using  TASAT 
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BRDFs  and  Eq.  (3).  Fig.  4  shows  simulated  brightnesses  for  this  satellite  during  the  fully-sunlit  portions  of  two 
LEO  terminator  passes  over  the  Maui  AMOS  site  and  the  South  American  SOAR  sites.  The  analyses  presented 
below  employ  simulations  from  these  two  ground-based  sites,  including  the  passes  shown  in  Fig.  4  and  other 
qualitatively- similar  additional  ones. 


Figure  3.  Schematic  illustrations  showing  two  views  of  the  simple,  non-convex  satellite  used  for  spectral  radiance  simulations  and  testing  of  the 
MAE  and  RDE  analysis  methods,  comprising  four  distinct  components  color-coded  as  follows:  “buscube”  (purple),  “bustube”  (red),  “wings” 
(green),  and  “struts”  (blue).  See  Table  1  for  details  on  the  specific  material  and  BRDFs  covering  each  of  these  components.  For  size  scale,  the 
sides  of  the  purple  cube  measure  1  m,  and  the  red  tube  1 .5  m  in  length. 


Component 
Index,  j 

Short 

Name 

Component 
color  in  Fig.  3 

BRDF 

Type 

TASAT 

BRDF 

Material  Description,  finish  and  other 
details 

1 

Bustube 

Red 

Maxwell  Beard 

0046 

Aluminum  Alloy,  2024-T3 

2 

Wings 

Green 

Maxwell  Beard 

0020 

Solar  Cell,  Silicon,  Sun  Side 

3 

Struts 

Blue 

Maxwell  Beard 

0029 

Aluminum  Alloy,  5456-H116,  Mill  Finish 

4 

Buscube 

Purple 

Maxwell  Beard 

0069 

Insulation,  Polyimide  Fiber,  Woven 

4.1  Results  for  Method  1:  Material  Abundance  Estimation  (MAE) 

The  MAE  analysis  method  has  previously  been  applied  to  a  convex  object  [1],  along  with  tests  of  truth  retrieval 
from  noise-free  scans,  as  well  as  discussions  of  the  limitations  of  noisy  measurement  data  and  BRDF  databases  that 
are  either  incomplete  and/or  possess  superfluous  materials  that  are  not  found  anywhere  on  the  observed  satellite. 
These  general  results  and  conclusions  apply  equally  well  to  the  non-convex  satellite  analyzed  here,  and  will  not  be 
restated  in  the  interest  of  brevity.  Instead,  this  analysis  focuses  on  results  from  simulations  of  noisy,  multi-site, 
ground-based  observations  of  the  satellite  in  LEO,  using  the  same  plot  formats  and  notation  introduced  in  [1]. 

4.2  Results  for  Method  2:  Reflectance  Distribution  Estimation  (RDE) 

The  RDE  method  has  a  much  grander  objective  than  the  previously-described  MAE  method,  in  that  it  seeks  to 
derive  the  complete  reflectance  distribution  functions  for  the  entire  set  of  satellite  components,  rather  than  just  a 
simple  mixture  of  materials  for  each.  Perhaps  the  most  important  lesson  learned  from  the  development  and  testing 
of  these  two  analysis  methods  is  that  the  RDE  method,  when  compared  to  the  MAE  method,  requires  a  much  larger 
amount  of  data  with  significantly  greater  geometric  diversity  to  converge  to  accurate  results.  Evidently,  to  derive 
accurate  specular  reflection  components,  the  observations  must  sample  each  satellite  component  glinting,  ideally 
from  several  perspectives.  In  other  words,  the  data  must  sufficiently  sample  the  “sharp  peaks”  of  the  BRDF 
surfaces,  such  as  those  shown  in  Fig  1.  Achieving  this  requires  a  tremendous  quantity  of  diversity  in  observational 
geometry  (i.e.,  richly  varied  viewing  and  illumination  perspectives). 

In  fact,  it  turns  out  that  the  RDE  method  can  suffer  numerical  instability  when  provided  with  input  data  of  limited 
geometrical  diversity.  The  reason  for  this  is  that  it  often  has  no  data  to  constrain  the  albedos  of  the  sharpest  specular 
components  in  the  retrieved  BRDF  expansion  series.  In  other  words,  the  amplitude  of  the  specular  BRDFs  shown 
on  the  lower-right  of  Fig.  2  can  grow  unphysically  large,  simply  because  few  (or  no)  observations  provided  data  that 
sampled  that  specular  peak.  This  manifests  itself  in  the  component  albedos,  aj9  approaching  unity,  because  the  most- 
specular  component  albedos,  a^m,  have  grown  to  large  values.  However,  there  is  a  numerical  means  of  addressing 
this  instability.  Specifically,  imposing  the  inequality  constraints  of  Eq.  (15)  employs  Lagrange  multipliers  [6,  7], 
one  for  each  BRDF  expansion  function  m.  These  multipliers  can  be  inspected  to  determine  if  this  kind  of  instability 
in  the  most  specular  BRDF(s)  has  manifested  itself  in  the  solution.  If  that  is  found  to  be  the  case,  then  these  affected 
BRDFs  can  be  eliminated  from  the  expansion  function  set,  the  corresponding  aJ>m  values  set  to  zero,  and  the  process 
repeated  until  a  numerically  stabilized  solution  is  found.  The  resulting  solution  set  {aj>m}  can  then  be  inserted  into 
Eqs.  (11)-(13)  to  calculate  the  best- fit  BRDF  for  each  satellite  component,  pj  (2,  n,  o,  s). 


Table  1.  Materials  and  BRDFs  used  for  satellite  components  in  the  spectral  radiance  simulations. 
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Figure  4.  Simulated  ground-based  observation  of  the  non-convex  satellite  shown  in  Fig.  3  flying  in  the  orbit  of  SSN  22176  making  passes  over 
SOAR  in  Chile  (left)  and  AMOS  on  Maui  (right)  on  2012  Jan  01.  The  simulation  assumes  nadir/velocity  stabilization,  with  the  red  tube  of  the 
satellite  pointed  along  the  nadir  direction  and  the  solar  panel  struts  aligned  along  the  velocity  direction.  The  dark  panels  at  the  top  shows 
temporal  signatures  along  with  renderings  of  the  cube  reflecting  sunlight  (upper  right),  its  position  along  the  ground -track  (upper  left)  and  the 
whole-body  I-band  (k  =  0.8  pm)  brightness  in  range-normalized  stellar  magnitudes  (plotted  for  fully  sunlit  conditions).  The  bottom  panels  show 
surface  plots  of  the  associated  multi -band  OCS  values  for  each  pass,  used  as  input  to  the  MAE  and  RDE  analyses  [1].  The  simulations  used  in 
this  analysis  employ  9  bands  (0.2,  0.4,  0.6,  0.8,  1.0,  and  1.2  pm)  and  assume  10%  Gaussian  noise  on  the  OCS  measurements. 

Fig.  5  shows  the  results  of  the  MAE  method  applied  to  progressively  increasing  amounts  of  ground-based  data.  The 
three  panels  show  the  results  from  one,  three,  and  six  terminator  passes  from  left  to  right.  It  demonstrates  that 
adding  more  data  to  the  process  leads  to  progressively  improved  MAE  results.  It  also  demonstrates  that  the 
smallest-sized  component  of  the  satellite  (the  struts  in  Fig.  3  and  Table  1),  which  naturally  reflects  less  light  than  the 
other  larger  components,  requires  the  most  data  in  order  to  converge  to  accurate  material  abundance  results. 


2  4  6  8  10 

Matenrt  Number 


Figure  5.  Estimated  material  abundances,  fj>m,  derived  by  combining  using  increasing  numbers  of  passes  of  data.  The  three  panels  show  the 
results  from  one  SOAR  ground-based  pass  (left),  three  passes  spanning  one  day  (1  SOAR  +  2  AMOS,  center),  and  six  passes  spanning  two  days 
(2  SOAR  +  4  AMOS,  right).  Each  plot  shows  the  four  satellite  components  (described  in  Table  1)  along  the  vertical  axis  and  twelve  candidate 
materials  along  the  horizontal  axis.  A  diagonal  pattern  red-squares  (fp„  =  1  along  the  diagonal  and^m  =  0  elsewhere)  in  the  first  four  columns  of 
each  plot  represents  perfect  truth  retrieval  [see  1].  As  can  be  seen  in  this  progression,  adding  more  data  to  the  process  leads  to  progressively 
improved  MAE  results  that  approach  perfect  truth  on  the  right  panel. 

Fig.  6  shows  the  best-fit  BRDFs  (along  with  known  truth),  derived  using  an  RDE  analysis  of  8  days  worth  of  data 
from  the  two  sites  (a  total  of  18  passes)  on  the  nadir-velocity  stabilized  model  satellite.  A  close  inspection  of  the 
BRDF  surfaces  and  albedo  spectra  plotted  in  Fig.  6  clearly  shows  that,  even  with  this  large  number  of  passes,  the 
method  really  has  not  converged  to  known-truth  reasonable  accuracy  (recall  that  the  MAE  method  converged  to 
reasonable  accuracy  using  just  6  passes).  Evidently,  even  eighteen  passes  observed  from  these  two  ground-based 
sites  still  does  not  provide  sufficient  geometric  diversity.  Two  ways  to  achieve  more  geometric  diversity  for  such  a 
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stabilized  satellite  would  be  to  add  more  ground-based  sites,  and/or  ever  larger  numbers  of  observed  passes. 
However,  this  research  indicates  that  a  far  more  efficient  way  is  to  employ  one  (or  more)  space-based  sensors,  which 
provide  more  dramatic  changes  in  viewing  and  illumination  perspectives  on  much  shorter  time -scales  than  can  be 
reasonably  achieved  even  with  a  very  large  number  of  ground-based  sensors.  Preliminary  analysis  indicates  that 
after  including  two  or  three  ground-based  sites  with  good  global  coverage  (especially  in  latitude),  adding  a  single 
space-based  sensor  provides  by  far  the  best  means  of  enriching  observational  diversity  when  observing  satellites 
already  in  orbit. 


Stabilized  Satellite 


=   •  ?  i 

Eighteen  Passes 
(7  SOAR  + 

11  AEOS) 

(8  days) 


Figure  6.  Estimated  BRDFs  derived  by  combining  using  18  passes  of  data  on  the  nadir-velocity  stabilized  satellite.  The  four  outer-most  panels 
show  the  known-truth  and  reconstructed  BRDF  surfaces  (top  left  and  right,  respectively)  and  the  known-truth  and  reconstructed  albedo  spectra 
(bottom),  plotted  in  black  curves  and  as  the  red  points,  respectively.  The  red  points  represent  the  best-estimate  total  albedos  for  each  of  the  9 
bands  in  the  simulated  data  set,  and  error  bars  account  for  the  “best-fit  uncertainty”  of  each  point,  but  not  for  any  deficiencies  in  the  BRDF  series- 
expansion  model.  The  pink  points  in  each  plot  show  the  estimated  Lambertian  fraction  of  the  albedo.  Each  panel  shows  the  results  for  the  four 
satellite  components  described  in  Table  1 ,  as  shown  by  the  arrows  emanating  from  the  central  rendering  of  the  model  satellite.  As  can  be  seen, 
even  a  total  of  eighteen  complete  terminator  passes,  compiled  over  eight  days  of  observations  from  the  two  ground-based  sites,  do  not  provide 
sufficient  geometric  diversity  to  converge  accurately  to  the  known-truth  BRDFs  for  this  nadir-velocity  stabilized  satellite. 

Changing  an  object’s  attitude  can  also  provide  increased  geometric  diversity.  Three-axis  satellite  stabilization,  such 
as  the  nadir-velocity  attitude  used  in  the  previously-analyzed  simulations,  tends  to  significantly  limit  geometric 
diversity  because  it  produces  very  slow  changes  in  the  relative  orientation  of  the  sensor(s)  and  the  Sun.  Generally, 
an  observer  would  not  be  able  to  change  an  orbiting  object’s  attitude  without,  at  the  very  least,  cooperation  from  the 
satellite’s  owner/operator.  However,  an  object’s  attitude  can  easily  be  changed  in  a  laboratory  observation 
environment,  or  in  simulation  studies.  Laboratory -based  tests  of  the  RDE  method  would  therefore  benefit  greatly 
by  using  multi-perspective  observations,  achieved  perhaps  by  mounting  the  object  on  an  articulating  robot  arm,  or 
on  a  rotating  platform.  Along  these  lines,  our  simulation  studies  indicate  that  rotational  motion  naturally  provides  a 
great  deal  of  geometric  diversity,  and  significantly  enhances  the  accuracy  of  RDE  method  analyses.  Furthermore, 
more  complicated  rotational  motion  produces  greater  diversity.  So  an  object  in  a  complex,  three-axis  rotation  state 
(e.g.,  precession+nutation)  would  generally  produce  more  diversity  than  an  object  in  a  stable  spin  (i.e.,  a  single-axis 
state). 


Fig.  7  shows  the  estimated  BRDFs  (along  with  known  truth)  derived  using  an  RDE  analysis  of  the  same  8  days 
worth  of  simulated  data  from  the  two  ground-based  sites  as  shown  in  Fig.  6,  but  with  the  satellite  in  a  three-axis 
rotation  state,  with  the  rotation  rates  about  each  axis  in  the  range  of  1  to  3  degrees  per  second,  so  that  the  object 
undergoes  several  complete  rotations  in  all  axes  during  each  observed  pass.  As  can  be  seen  in  Fig.  7,  this  dramatic 
increase  in  geometric  diversity  allows  the  RDE  method  to  converge  close  to  known-truth  with  reasonable  (but  not 
perfect)  accuracy.  The  remaining  differences  between  the  known-truth  and  best-fit  BRDF  surfaces  (top  plots  in  the 
four  outer-most  panels  of  Fig.  7),  as  well  as  in  the  known-truth  and  best-fit  albedos  (the  black  curves  and  red  points 
in  the  bottom  plots  of  the  four  panels)  are  likely  due  to  a  combination  two  inadequacies  of  the  current  RDE 
implementation:  1)  that  too  few  terms  were  used  in  the  BRDF  series  expansion,  and/or  2)  that  the  Lambertian  + 
Cook-Torrance  analytical  functions  themselves  cannot  capture  all  of  the  relevant  aspects  of  real-world  BRDFs. 
Resolving  these  remaining  discrepancies,  and  determining  more  specifically  what  constitutes  sufficient  geometric 
diversity  for  the  CDE  method  to  converge  accurately,  requires  further  simulation  and/or  laboratory  studies. 
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Figure  7.  Estimated  BRDFs  derived  from  data  on  a  rotating  satellite,  shown  in  a  format  similar  to  that  of  Fig.  6.  The  added  geometric  diversity 
naturally  provided  by  the  object’s  rotational  motion  allows  the  RDE  analysis  to  converge  much  more  accurately  to  the  known-truth  BRDFs.  The 
remaining  differences  between  the  known-truth  and  best-fit  BRDF  surfaces  (top  plots  in  the  four  outer-most  panels),  as  well  as  in  the  known-truth 
and  best-fit  albedos  (the  black  curves  and  red  points  in  the  bottom  plots  of  the  four  panels)  are  due  to  a  combination  of  inadequacies  in  the 
number  of  terms  used  in  the  BRDF  series  expansion,  and/or  the  in  the  quality  of  the  BRDF  basis  functions  themselves. 


5  CONCLUSIONS  AND  FUTURE  WORK 

This  paper  presents  the  theory  required  to  retrieve  satellite  surface  properties  from  temporal  sequences  of  whole  - 
body,  multi -band  brightness  measurements,  focusing  on  two  new  analysis  methods.  The  first,  material  abundance 
estimation  (MAE),  determines  the  abundances  of  materials  covering  the  components  of  the  observed  satellite.  The 
second,  reflectance  distribution  estimation  (RDE),  determines  the  BRDFs  of  the  satellite’s  components.  The  MAE 
method  can  suffer  from  the  limitations  of  the  BRDF  database  of  candidate  materials  that  it  requires,  making  it 
inappropriate  for  unknown  or  aging  satellites.  The  RDE  method  was  developed  specifically  not  to  need  such  a 
database,  but  instead  estimate  BRDFs  for  each  of  the  satellite’s  components  using  a  series  expansion  approach.  The 
RDE  method  requires  data  with  significant  geometric  observational  diversity,  in  order  to  converge  with  reasonably 
accuracy.  Employing  multiple  ground-based  sites  can  help  provide  such  diversity,  but  a  much  more  efficient 
approach  would  be  to  use  even  a  single  space-based  sensor.  Further  analysis  needs  to  be  conducted  to  determine  if 
achieving  sufficient  diversity  for  stabilized  satellites  is  even  possible,  and  if  so,  if  it  requires  a  prohibitive  number  of 
ground-  and/or  space-based  observations. 
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